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ABSTRACT 

With more and more extrasolar planets discovered in and around binary star systems, ques- 
tions concerning tlie determination of the classical Habitable Zone arise. Do the radiative and 
gravitational perturbations of the second star influence the extent of the Habitable Zone sig- 
nificantly, or is it sufficient to consider the host-star only? In this article we investigate the 
implications of stellar companions with different spectral types on the insolation a terrestrial 
planet receives orbiting a Sun-like primary. We present time independent analytical estimates 
and compare these to insolation statistics gained via high precision numerical orbit calculations. 
Results suggest a strong dependence of permanent habitability on the binary's eccentricity, as 
well as a possible extension of Habitable Zones towards the secondary in close binary systems. 

Subject headings: Habitable Zone — binary stars — S-Type 



Introduction 



Fueled by the successes of cur rent transit 



observation incentives like KEPLER (Welsh et al 



2012: 'Bo rucki k Kochl201lh and CoRoT (jBaglin et al 
2009i : iTinglev et al.l 120111) the quest for discover- 
ing the first Earth-twin has lead to a consider- 
able cross-disciplinary interest in the interplay 
between stellar and planetary properties to pro- 
duce habitable worlds. Even-though opinions 
differ on what exactly to look for in a system 
harboring a terrestrial planet in order to de- 



clare it 'habitable' ( see e.g . iBuccino et al.l (2006), 



Kaltenegger et al.| |2007l) . ISelsis et al.l 



20071) . 



Lammer et al.l (120091) the c l assica l assumption 



investigated bv iKasting et al.l (|1993^ . i.e. the ca 



pacity for water to stay liquid on the planet's 
surface, may still be considered a prerequisite for 
the development and sustainability of comple x 
life as we know it ([Kaltenegger fc Sasselovll2Qlir) . 
This entails restrictions on planetological charac- 
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teristics, such as mass, atmospheric and bulk com- 
position, and sets limits to the ho st star's activ- 



ity as well as radiation properties (|Lammer et al 



20091 ). Dynamical considerations are of equal im 
portance, since changes in orbital stability, or ex- 
treme variations in insolation due to large plan- 
etary eccentricities (er, > 0.7) may also result i n 
a hostile environment ( Williams fc PollardI 2002 ) . 
It is therefore only natural that one would look 
for a copy of our Solar System, when searching for 
habitable worlds. Yet, the study of exoplanetary 
systems so far has clearly shown that a broader 
perspective is required. 

In fact, up to 70% of all stellar systems in our 
galaxy may not be single but multi-stell a r sys - 
tems (e.g. iKiseleva-Eggleton fc EggletonI (|200ll ) 
and references therein). Together with the ap- 
proximately 60 planets that have alread y been dis- 
cover ed in systems harboring two stars (j Schneider 
20111 ) this suggests that binary and multiple star 
systems should not be ignored in the search for 
habitable worlds. 

Investigations of environments that permit 
planetary formation in binary star systems have 
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progresse d rapid ly over the last decade (see e.g. 
Thebauld (j201l[ ) and references therein). Even- 



though important questions regarding the early 
phases of planet formation in binary star sys- 
tems - especially the transition from planetesi- 
mal to planetary embryos - still remain to be 
answered, late stage formation scenarios for ter- 
restrial planets in such environmen ts are available 
(iWhit mire et al.l 1998; Haghighipo ur fc Ravmondl 
[2OO7I: iHaghighipour et al.l 120101 ). Since previous 
studies did consider the extent of the classical 
Habitable Zone (KHZ) to be purely a function of 
the primary star 's luminosity and sp ectral type 
as introduced by iKasting et al. I (Il993h (hereafter 
KWR93), we aim to refine this definition to en- 
compass the gravitational and radiative influence 
of a second star. 

This article is structured as follows: After a 
short recapitulation of the main radiative aspects 
of habitability as defined in KWR93 section [5] in- 
troduces three exemplary binary-planet configura- 
tions, which will serve as test-cases for habitabil- 
ity considerations. Section [3] briefiy mentions the 
dynamical requirements which binary-planet con- 
figurations have to fulfill in order to ensure system 
stability. In section 2] the maximum radiative in- 
fluence of the second star on a terrestrial planet 
in the primary's HZ is estimated and compared 
to actual insolation simulations. The occurring 
differences are being investigated in the following 
section. Finally, generalized, analytical estimates 
are developed and compared to numerical simula- 
tions in sections [6] & [Zl and the results concerning 
the behavior of HZs in binary star systems are 
presented in section El A discussion of the results 
concludes this article. 

2. Binary-Planet Configurations 

Apart from planetological and dynamical as- 
pects, the insolation a terrestrial planet receives 
from its host star is naturally the main driver de- 
termining the extent of the HZ. When considering 
planets within a binary star system it is there- 
fore important to track the combined radiation of 
both stars arriving at the planet. KWR93 showed 
that not only the sheer amount of insolation, but 
also the spectral distribution is essential to es- 
timate limiting values for atmospheric collapse. 
In order to model the impact of diverse stellar 



spectral classes on an Earth- like planet's atmo- 
sphere, KWR93 introduced so-called effective ra- 
diation values. These measure the relative im- 
pact a comparable amount of radiation (e.g. 1360 
[Vl^/m^]) with different spectral properties has on 
a planet's atmosphere. Taking the effects of differ- 
ent stellar spectra into account is especially impor- 
tant in binary star systems with different stellar 
components. Tab. [1] reproduces the effective ra- 
diation values for the inner (runaway greenhouse) 
and outer (maximum greenhouse) edge of a single 
star's HZ as given in KWR93. Notice how the on- 
set of runaway greenhouse effects requires almost 
twice as much radiation for a spectral distribu- 
tion akin to FO class star s compared to MO spec- 
tral ty pes. Even though iKaltenegger fc Sasselov 
(j201l[ ) assume similar effective radiation values for 
M and K spectral classes, actual calculations have 
only been done for FO, G2 and MO ZAMS stars. 
For this reason, we first investigate the following 
three stellar configurations: 

i) G2 - MO ^ ~ 0.3 

ii) G2 - G2 /i ~ 0.5 

iii) G2 - FO /i ~ 0.6 

All stellar components are considered to be ZAMS 
stars, and /i = 1712/(1711 + 1712) denotes the bi- 
nary's mass ratio. A terrestrial planet is orbit- 
ing the Sun-like G2 host-star, hereafter referred 
to as primary. Such binary-planet configurations 
are considered to be of satellite type ( S-Type), 
i.e. the planet revolves around one star ( DvorakI 



19841 ). see Fig.[TJ In fact, most of the planets in bi- 
nary systems hav e been discovered to be of S-Type 
Schneideil 2011), e .g. the system Gamma Cephei 



Hatzes et al.l 2003 ). In the following, we choose 



binary systems with semi-major axes between 10 
and 50 AU for the comparison of numerical and 
analytical estimates on the extent of the HZ, as 
for closer binaries the HZs in G2-G2 and G2-F0 
configurations are considerably reduced due to dy- 
namical instability, whereas the qualitative differ- 
ences for results beyond 50 AU are small. How- 
ever, the methods presented in section[6]are viable 
beyond those limits, as long as the assumptions 
given in iGeorgakarakos (2003) remain valid. 

The terrestrial planet is assu med to be cloudless 
(jKaltenegger fc Sasselovl 120111 ) and orbits the G2 
star only, resulting in configurations i) and ii) to 
be classified as S-Type I (^ < 0.5) and in) as S- 
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Table 1 

Limiting radiation values for the inner (A) and outer {B) border of the HZ respectivel y 



IN UNITS of Solar constants (1360 [PF/m^]). The values were taken from IKasting et al. 



([19931) assuming a runaway greenhouse scenario for the inner limit, and a maximum 

greenhouse effect for the outer limit. 



Spectral Type 


A 


B 


FO 


1.90 


0.46 


G2 


1.41 


0.36 


MO 


1.05 


0.27 



Type II (^ > 0.5) respectively (see Fig.[T]). 

3. Dynamical Stability 

Binary-planet configurations that lead to highly 
chaotic orbits and eventual ejection of the planet 
cannot be considered habitable. Therefore dynam- 
ical stability of the system is a basic requirement 
for our considerations. In order to assess the dy- 
namical stability of a test planet's orbit within the 
given binary star systems, we applied results of nu- 
merical stability s tudies by Rabl fc Dvorak ( 19881) 
(hereafter RD88),'Pilat-LohinKer & Dvorak' (200^ 
(hereafter PLD02) and H olman fc Wie gcrt (199_^ 
(hereafter IIW99) determining stable zones in a 
planar, normalized system (binary semi-major 
axis ah — 1). In contrast to RD88 and IIW99 
who classified unstable orbits via ejections of test 
planets from the system, PLD02 applied the Fast 
Lyapunov Indicator (F LI) chaos detection method 
( Froeschle et al. 1997t) im plemented in a Bulirsch - 
Stoer extrapolation code (jBulirsch fc Stoeiill964[) . 
In Tab. [5] maximum planetary semi-major axes 
which still allow for a stable configuration in the 
normalized setup are shown for different binary 
orbits and mass ratios. Results gained via FLI 
by PLD02 are compared to those published in 
HW99. Even-though the critical values are very 
similar, the onset of dynamical chaos (PLD02) 
in G2-G2 and G2-F0 configurations appears for 
smaller planetary semi-major axes than predicted 
by the ejection criterion in IIW99. For the fol- 
lowing investigations the more conservative FLI 
stability estimates by PLD02 are used. 

4. Influence of the Secondary 

The main question in dealing with habitabil- 
ity in binary star systems is doubtlessly: "How 



does the second star affect the HZ around the 
primary?" Let us focus on the dynamical aspects 
first, as stability is a prerequisite to habitability. 
The results from the previous section dictate that 
not all combinations of a binary's semi-major axis 
ttb and eccentricity are viable, if orbital stability 
of planets in the primary's HZ is required. Fig. [31 
top left shows quadratic fits of the FLI stability 
data presented in Tab. [51 The secondary's allowed 
eccentricities are higher for planets orbiting the 
primary near the inner border of the KHZ, than for 
cases where the planet's semi-major axis is close to 
the KHZ's outer rim. Such dynamical restrictions 
set limits on how close the secondary can approach 
the planet, which will in turn limit its insolation 
on the planet. Applying analytical approxima- 
tions introduced later in this section Fig. [3l top 
right demonstrates that the smallest possible sep- 
aration between planet and secondary is always 
larger than the planet's aphelion distance to its 
primary. This lessens the second star's poten- 
tial radiative contribution considerably for all but 
close S-Type II configurations. In the latter case 
the enhanced luminosity of the secondary compen- 
sates larger distances to the planet. Using the 
minimum distances presented in Fig. [3l top right, 
one can estimate, that in a G2-G2 S-Type I system 
with Ob = 2QAU , the secondary's radiative infiu- 
ence on a planet started at the inner edge of the 
KHZ is in the order of only 10% of the primary's 
contribution. As a consequence, the secondary 
was often considered not to have a significant ra- 
diative impact on the extent of the primary's HZ 
(Whitm ire et al.l 1998 ': Haghi ghipour fc Ravmond 
[2007. : .Haghighipour et al...20l6[ l 

Such a so-called "single-star approach", how- 
ever, stands in contrast to numerical experiments 
presented in Fig. [H left. The effective insola- 
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Table 2 

Critical, planetary sem i-major axis in a no r maliz ed binary star system ai, — 1 with 

MASS- RATIO M AS STATED IN IHOLMAN &: WlEGERll (j 19991 ) COMPARED TO THE VALUES FOUND BY 
(jPlLAT-LOHINGER fc DvORAkI|2002| ') VIA FAST LyAPUNOV CHAOS INDICATOR (FLI). THE HIGHER THE 
BINARY'S ECCENTRICITY Cf,, THE SMALLER THE TEST-PLANET'S SEMI-MAJOR AXIS HAS TO BE TO PERMIT 
DYNAMICAL STABILITY. EVEN-THOUGH THE RESULTS ARE VERY SIMILAR, THE FLI LIMITS TEND TO BE 
MORE CONSERVATIVE FOR HIGHER MASS RATIOS AND BINARY ECCENTRICITIES. 







= 0.3 




= 0.5 




= 0.6 




FLI 


HW99 


FLI 


HW99 


FLI 


HW99 


0.0 


0.37 


0.37 


0.27 


0.26 


0.23 


0.23 


0.1 


0.29 


0.30 


0.25 


0.24 


0.21 


0.20 


0.2 


0.25 


0.25 


0.19 


0.20 


0.18 


0.18 


0.3 


0.21 


0.21 


0.16 


0.18 


0.15 


0.16 


0.4 


0.18 


0.18 


0.15 


0.15 


0.12 


0.13 


0.5 


0.13 


0.14 


0.12 


0.12 


0.09 


0.10 


0.6 


0.09 


0.11 


0.08 


0.09 


0.07 


0.08 


0.7 


0.07 


0.07 


0.05 


0.06 


0.05 


0.05 


0.8 


0.04 


0.04 


0.03 


0.04 


0.02 


0.035 


0.9 


0.01 




0.01 




0.01 





tion curves shown were generated solving the full 
Newtonian Three-B ody Problem (3BP) numeri- 
cally via Lie- Series ( Hanslmeier fc Dvor ak _ 1984t 



Eggl fc Dvora3l2nin[ ) and Gauss-Radau teverharti 



19741 ) integrators, where the actual amount of ra- 



diation the planet receives was calculated for each 
integration step. One can see that a planet started 
at the inner edge of the KHZ in a G2-G2 S-Type 
binary experiences an increase of more than 30% 
in insolation compared to a planet on a circular or- 
bit with corresponding semi-major axis around a 
single G2 star. Consequently, an Earth-like planet 
would spend a considerable time outside the clas- 
sical, circular HZ in a G2-G2 configuration. Is, 
therefore, the secondary's radiative impact more 
important than assumed? The fact that the vari- 
ations in insolation perfectly correlate with the 
dynamical evolution of the planet's eccentricity 
(see Fig. 21 right) permits an alternative expla- 
nation: changes in the planetary orbit induced 
via gravitational perturbation by the second star 
might be responsible for the increased insolation 
values. Even-though planetary and binary orbits' 
semi-major axes are n ot expected to sh ow any 



secular variation (e.g. Harrington! ( 19681 )). it is 



known that even a distant companion would in- 
ject some eccentricity into the orbit of two bod- 
ies rev olving around their comm on center of mass 
(e.g. iMazeh fc ShahamI ^7^; \g eor gakarakosi 
(|2n02f )'l. Elevated planetary eccentricities would 
entail smaller periastron distances - allowing for 



increased insolation by the primary. 

In order to draw a clearer picture on whether 
the perturbation induced rise in the planet's eccen- 
tricity or the secondary star's radiative influence 
are the dominant factors causing the increased in- 
solation onto the planet encountered in Fig.|4l left, 
we will make three assumptions which our analyt- 
ical estimates will be based on: 

a) the binary-planet system is coplanar 

b) stellar luminosities L are constant on timescales 
of the system's secular dynamics. 

c) stellar occultation effects are negligible. 

Since we assume coplanar orbits, we wi ll make ex- 
tensiv e use of the analytic expressions in Georgakarakod 
(|200l I2OO5) to calculate the time averaged 
squared planetary eccentricity (ep)^. In Appendix 
IB] we will derive estimates on the planet's max- 
imum injected eccentricity e^° ^. Unlike the re- 



cen t ansatz by iGiuppone et al.l (|2011l ) and ear- 
lier Thebault et al.l (|2006 ). where the eccentricity 
evolution of a planet in a st ellar binary was mod- 
eled by empirical formulae, Georgakarakod ( 2003 , 
20051) derived an entirely analytical formula for 
the eccentricity of the inner - in this case the 
planet's - orbit of a hierarchical triple system, 
which is valid for a wide range of mass ratios and 
initial conditions. Together with the estimates 
in Appendix [B] the formulae presented 
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in Georgakarakos ( 2003 . 2005f) are an analytical 
extension to the the first order with respect to 
the perturbing mass secular perturbation theory, 
as given for example in the longstanding work of 
Heppenheiineil (|l978[) . 

Assumptions b) and c) are reasonable for well 
separated binary stars, where both components 
are on the main sequence. The combination of 
the dynamical stability results presented in Fig. [31 
top left, and the analytic expressions for the in- 
jected planetary eccentricity allow us to estimate 
not only the minimum distances between the sec- 
ondary and the planet as shown in Fig. [H top 
right, but also the maximum contributions to the 
planetary insolation from the primary and the sec- 
ondary star respectively. 

Fig. El bottom left shows the primary's effective 
insolation onto a terrestrial planet during its clos- 
est perihelion passage d = ap{l — e™""^). Here, the 
largest dynamically possible perturbation by the 
secondary is considered. In order to separate the 
secondary's radiative and gravitational effects, the 
second star's radiation has been excluded in this 
plot. Given a binary semi-major axis of lOAU, 
Earth-analogues started at the inner edge of the 
KHZ permit considerably higher binary eccentric- 
ities (Fig.|3l top left) compared to planets near the 
outer edge. Therefore, higher insolations for plan- 
ets near the inner edge of the KHZ are expected, 
since the injected planetary eccentricities are cou- 
pled strongly to the binary's eccentricity, see Ap- 
pendix [Bl This effect abates around ai, = 20AU, 
where the primary's insolation becomes almost 
equal for planets started near the inner and outer 
rims of the KHZ. For binary semi- major axes be- 
yond Ub = 20AU planetary orbits closer to the 



I star 2 



Starl 



Fig. 1. — Two examples of S-Type motion, i.e. the 
plane t orbits only one binary-component ([Dvorak 
19841 ): left: S-Type I {p < 0.5), the more mas- 
sive star is the planet's host (primary), right: 
S-Type II {fi > 0.5), the less massive star is the 
planet's host. 



secondary are more severely perturbed, than the 
ones near the inner edge of KHZ, as the binaries' 
maximum allowed eccentricities become similar for 
inner and outer borders KHZ and the perturbation 
is merely dependent on the binary to planet period 
ratios. 

For planets on eccentric orbits in a G2-G2 con- 
figuration with Qb = 20AU, the primary's insola- 
tion can increase up to almost 50% compared to 
circular orbits sharing the same semi-major axis. 
In contrast, even with the planet being as close to 
the second star as dynamically possible, the sec- 
ondary's maximum radiative input only accounts 
for an additional 10% in the same setup. This 
can be seen in Fig. [31 bottom right. Hence, the 
primary is the main source of the additional in- 
solation in S-Type I systems, which can only be 
explained via changes in the planet's eccentricity. 
Solely for S-Type II systems like G2-F0 binaries, 
the secondary's radiative contributions to plane- 
tary insolation can become comparable to the pri- 
mary's (cf. Fig.[3l bottom). 

5. Consequences for the Habitable Zone 



In the previous section, the rise in the planet's 
eccentricity was identified as a main driver for the 
strong variance in simulated insolation curves in 
S-Type I systems. Eccentric planetary orbits en- 
tail strong variations i n insolation over an orbita l 
period. Nevertheless, [Williams fc PoUardl (|2002h 
concluded, that this might not necessarily be pro- 
hibitive for habitability. As long as the average 
insolation values lie within habitable limits and 
Cp < 0.7, the atmosphere should be able to act 
as a buffer, preventing immediate carbon freeze- 
out or instantaneous water evaporation. In order 
to distinguish between cases, where the planet re- 
mains within the radiative boundaries of the KHZ 
for all times, configurations where the planet is 
"on average" within the HZ, and non habitability, 
we introduce the following three categories: 



Permanently Habitable Zone (PHZ): The 

PHZ is the region where a planet always 
stays within the insolation limits {A, B) of 
the corresponding KHZ, i.e. 

A > Seff > B (1) 

where A denotes the inner, and B the outer 
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effective radiation limit for a given spectral 
type, see Tab. [T] 

Extended Habitable Zone (EHZ): In con- 
trast to the PHZ - where the planet stays 
within the KHZ for all times - parts of the 
planetary orbit lie outside the HZ due to e.g. 
the planet's eccentricity. Yet, the binary- 
planet configuration is still considered to be 
habitable when most of its orbit remains 
inside the HZ boundaries: 

{Seff),+<7<A A {S,ff)^-a>B (2) 

where {Seff)f denotes the time-averaged ef- 
fective insolation from both stars and the 
effective insolation variance. 

Averaged Habitabl e Zone (AHZ): F o llowin g 
the argument of IWilliams fc PoUardl (|2002l ). 
this category encompasses all configurations 
which allow for the planet's time-averaged 
effective insolation to be within the limits of 
the KHZ. 



A > {S^ff), > B 



(3) 



6. Analytical Estimates 



We now propose analytical estimates to achieve 
a classification of planetary habitability as was 
suggested in the previous section. The aim is to 
circumvent time consuming numerical integrations 
when global parameter scans are required to check 
systems for possible habitability. Even-though the 
following analytical estimates are presented utiliz- 
ing the Seff values developed in KWR93, more 
advanced atmospheric models for exoplanets can 
easily be introduced by exchanging the effective 
insolation values A and B. 

6.1. Estimates for the PHZ 

Let the second star move on a fixed Keplerian 
orbit with semi-major axes Ub and eccentricity et- 
Accordingly, the planet's orbit has the semi-major 
axis Gp and acquires a maximum eccentricity of 
gmax j^^g i^j^g secondary's gravitational pertur- 
bations (cf. Fig m right). This permits us to 
estimate the maximum and minimum insolation 
conditions for the planet to permanently remain 
within the KHZ: 



insolation minimum condition: Both, planet 
and secondary are assumed to be in apocen- 
ter position and opposition. The additional 
normalization of the stellar luminosities per 
solid angle (L — Lboi/{4:n)) with regard to 
the respective outer insolation limits for each 
star (Bi, B2) ensures that different spectral 
properties are taken into account. 



1 < 



|^(ap(l + e™--)) + 
||(ab(l + e5) + ap(l + e™-))- 



insolation maximum condition: Again the lu- 
minosities are normalized, but this time with 
regard to the inner insolation limits (Ai, 
A2). Since we consider S-Type I and II sys- 
tems, it is possible that the secondary at 
pericenter may produce a higher insolation 
on the planet than the primary star. If this 
is the case, the maximum insolation configu- 
ration will have the planet at apocenter with 
regard to the primary, and the secondary at 
pericenter. 



^(ap(l-e 



max 
P 



1 > max < 



^(ap(l-|-e^--)) ' 
[ +^(a,(l-e,)-ap(l + e™-))" 



6.2. Estimates for the AHZ 

The combined stellar insolation Stot on the 
planet is, of course, a function of time. In order to 
calculate time averaged insolation values, we will 
use that: 



{Stot)t = (Si), + (52), 



(4) 



where ("S*!), is the time average of the planetary 
insolation due to its host star, and (S'2)( the time 
averaged contribution of the second star. Let us 
focus on the two-body problem planet - host star 
first: Here, {Si)^ = {Li/5'^{t))^ where 5{t) denotes 
the scalar distance between planet and primary. 
The average insolation a planet on an unperturbed 
Keplerian orbit experiences can be calculated us- 
ing the planet's angular momentum h = 5^ f , f 
being the true anomaly. 



{Si)t 



P 



dt 



(5) 
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Lin 

~r 

This expression states that the time-averaged in- 
solation a planet receives over one orbit depends 
only on the star's luminosity, the planet's mean 
mot ion n and its orbital angular momentum (see 
e.g. ISeager ( 2010[ ). p. 18). Now let us construct a 
circular orbit so that hdrcie = heiupse- A planet 
moving on such an orbit with 'equivalence radius' 
r = ap{l — {elp)t) will experience the same amount 
of insolation per unit time as a planet on an el- 
liptic orbit sharing the same angular momentum. 
The advantage of considering "equivalent circu- 
lar orbits" is that the insolation remains constant 
for any orbital position of the planet. Of course, 
the reduction of elliptic orbits to circular orbits 
with common angular momentum will decrease 
the planet's orbital period, as r < a and r ~ a 
only if (ep)( = 0. However, as we have chosen the 
'equivalent circular orbit' to have the same angu- 
lar momentum, it will share the same constant rate 
of change of insolation with the true orbit. There- 
fore, even an average over the longer period of the 
elliptic orbit will still yield valid results. 

We will now apply the same line of argumenta- 
tion to construct the equivalence radius R for the 
secondary. This allows us to extend the pseudo- 
static radiation environment to include the aver- 
age radiative influence of the secondary. Since the 
distinction between averaged and initial eccentric- 
ity is small for the secondary - its secular variance 
is negligible for the cases investigated - we can 
safely use the secondary's initial eccentricity to 
calculate R. The suggested configuration is shown 
in Fig. [21 The secondary circles the primary in a 
fictitious orbit with equivalence radius R, and the 
planet orbits the primary in a circle with radius r. 

The insolation on the planet depends on the rel- 
ative distances of the planet to both stars. In order 
to estimate the insolation on the planet caused by 
the secondary (S'2), we simply apply law of cosine, 
and average over all possible geometric configura- 
tions (Fig.E]): 

= i?2 + - 2ri?cos((/)) 
L2 



{S2) 



(6) 




Fig. 2. — 'Equivalence Orbits' of the planet (dot- 
ted) and the secondary (dashed) around the pri- 
mary. R denotes the equivalence radius of the sec- 
ondary's circular orbit sharing the same angular 
momentum with its actual elliptic orbit (continu- 
ous). The planet's equivalence radius (r), the dis- 
tance between the time averaged secondary's and 
planet's positions (p) as well as the angle ((j)) op- 
posite to (p) are highlighted. In order to calculate 
the secondary's mean insolation on the planet, an 
averaging over 1/ p^{(j)) is required. 



L2 r 1 

27r io R^ + r"^ - 2rR cos(0) 
L2 



(i?2 _ r2) 



if R> r 



where p is the planet's distance to the secondary, 
and (j) denotes the angle between the distance vec- 
tors of planet and secondary to the host star. Since 
we do not allow for orbit crossings of the planet 
and the secondary, R > r will always hold. Con- 
sequently, the total, time averaged insolation onto 
the planet is given by: 



{Stot)t 



Li 

,2 



r2 (i?2_^2) 



(7) 



Equation (O does not yet take the different spec- 
tral properties of the binary's components into ac- 
count. Therefore we note the following conditions 
for the planet's averaged insolation being within 
habitable limits: 

average insolation minimum condition: 

L2 1 



1< 



Bi r2 



B2 i?2 - r2 
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average insolation maximum condition: 



Li 1 



1 



based on Gauss R adau quadrature (lEverhart 
1974) and Lie Series (iHanslmeier fc Dvoraklll984 



Eggl fc Dvorak 2010i) to determine the actual po- 



Here, the indices 1, 2 indicate the respective star's 
KHZ boundary values A and B, R = ab{l — 
el) and r = ap{l — (Cp)*) where the averaged 
square d planetary eccentrici ty was calculated fol- 
lowing Georgakarakod ( 20051 ). 



6.3. Estimates for the EHZ 

Having derived the insolation time averages in 
the previous section, the expected insolation vari- 
ance (cr^) still remains to be determined. 



= {Sfot)t - {Stotft 



(8) 



Using equation ([7]), and analytic estimates for 
i^totlt: the effective insolation variance can be cal- 
culated as follows (see Appendix 



4 



L 



^,lj(-l + 3(e^)-3(e^}2 + (e2)3) 

+ Xfr^y ^"^P'y 2 2 
2L1L2 



XiX2(r4-r2i?2) 



where Xi g {Ai,Bi} and the index i denotes the 
respective star. The minimum and maximum con- 
ditions for a planet to be within the EHZ are given 
as: 

extended insolation minimum condition: 



1 < {Seff,B)t - CTB 

extended insolation maximum condition 



1 > {Seff,A)t + f^A 

7. Reliability of Analytical Estimates 



sitions of both stars with respect to the Earth-like 
planet. Assuming that the stellar luminosities 
will not change significantly on the timescale of 
the planet's secular period in eccentricity, good 
approximations of insolation patterns can be ob- 
tained in such a way. Fig. [6] shows a comparison 
of analytic habitability classifications versus re- 
sults gained via numerical orbit integration and 
direct insolation calculation. The setup consists 
of a terrestrial planet (1 M®) in S-Type orbit 
around a G2 host star, with three different spec- 
tral types as secondary: FO (top), G2 (mid) and 
MO {bottom). The terrestrial planet was started 
on circular orbits with semi-major axes between 
O.QAU < Up < 2AU with the secondaries' semi- 
major axes being Ob = 50 AU. The time span of 
the numerical integrations encompassed at least 
two secular periods in the planet's eccentricity. 
It is evident that for small binary eccentricities, 
all three types of HZ coincide well with the bor- 
ders defined in KWR93 indicated by the vertical 
lines at 0.84 and 1.67 AU. For > 0.1 how- 
ever, a splitting into the HZ categories defined 

in sections [5] and [6] becomes eminent. The PHZ 

( 1 — A /l — (e2) (1 + (e^)) ] (black - blue online) shrinks considerably with 

growing eccentricity of the binary's orbit. This is 
due to the perturbation induced elevation of the 
planet's eccentricity. In contrast, the region de- 
fined in KWR93 is best approximated by the AHZ 
(light grey - yellow online), which remains virtu- 
ally unaffected by the secondary's eccentricity. 
In this setup all analytically calculated HZs are 
in excellent agreement with the numerical ones. 
Only close to the stability limit (shaded region - 
purple in the online version) the correspondence 
between simulation and analytical estimates de- 
creases. This can be seen more clearly when the 
secondary's influence becomes stronger, e.g. in 
the cases of at = 10 AU, see Fig. [T] In general the 
analytical approach is producing more conserva- 
tive results compared to the numerical data (cf. 
Figs. [5] & [71). In order to determine whether these 
deviations in the determination of the PHZ are 
due to 



(9) 



In order to test the reliability of the analyt- 
ical estimates presented in section |51 we used 
high precision numerical integration methods 



a) inaccurate analytical estimates for e^"^, 

b) insufficient time resolution in numerical sim- 
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Table 3 

Percentages of planetary orbits classified identically via a numerical simulations and 
analytical estimates as presented in section [sj three binary component configurations 
have been investigated, the reference classifications were extracted from numerical 
orbit integrations and insolation simulations. 



[AU] 




G2-M0 [%] 






G2-G2 [%] 






G2-F0 


[%] 




a.b 


Total 


PHZ 


EHZ 


AHZ 


Total 


PHZ 


EHZ 


AHZ 


Total 


PHZ 


EHZ 


AHZ 


10 


95.9 


97.4 


99.8 


98.3 


94.4 


97.4 


99.2 


98.0 


93.6 


95.8 


98.5 


98.2 


20 


98.8 


99.3 


99.5 


99.8 


98.5 


99.2 


99.6 


99.6 


98.5 


99.5 


99.4 


99.6 


30 


99.0 


99.5 


99.7 


99.9 


99.2 


99.7 


99.6 


99.8 


98.9 


99.8 


99.4 


100.0 


40 


99.2 


99.5 


99.6 


99.9 


99.3 


99.9 


99.6 


100.0 


99.0 


99.8 


99.5 


99.8 


50 


99.2 


99.6 


99.7 


99.9 


99.4 


99.7 


99.7 


99.9 


99.4 


99.8 


99.7 


99.9 



ulations with regard to determining e™"^ 
values, or 

c) insufficient total integration time to reach 
minimum and maximum insolation condi- 
tions, 

we constructed semi-analytical PHZs using numer- 
ically determined e™""^ values in the analytic equa- 
tions to determine the PHZ presented in section 
16.11 The borders of the semi-analytically derived 
PHZs are depicted as white dashed-dotted lines 
in Figs. [5] & [71 As they are nearly identical with 
the fully analytic estimates, we can exclude a) and 
b) , which would have lead to significantly different 
results for semi-analytic and analytic approaches. 
Therefore, c) seems most likely to cause the differ- 
ences between the numerical and analytical PHZs, 
since encountering an exact alignment of planetary 
aphelion and secondary perihelion at the moment 
where tp = e™°^ may take far longer than two 
secular periods. As the computational efforts re- 
quired to ensure that a simulation's time resolu- 
tion as well as total integration time are sufficient 
to identify the correct PHZ boundaries are enor- 
mous, the necessity to have analytical methods at 
hand becomes evident. 

As far as EHZ and AHZ regions are concerned, 
clear differences can be seen in high perturbation 
environments close to the transition to instabil- 
ity (Figs. El & III shaded regions). In these cases, 
the authors favor the numerical results, as single 
configurations are not critical for the more statis- 
tically oriented measures. 

A quantitative overview of the correspondence 
between numerical and analytical results for all 
system configurations investigated is given in 



Tab. [21 Here, similar maps as presented in 
Figs, m & [71 were generated with a resolution of 
Aop = 0.01 AU, and Ae;, = 0.01 and evaluated 
statistically. In spite of their shortcomings in de- 
termining the PHZ, the numerical classification 
results have been used as reference values and 
are compared to the analytical estimates given in 
section [6l The total correspondence percentages 
are calculated as the number of all orbits below 
the stability limit that were classified identically 
via numerics and analyitcs, divided by the total 
number of orbits simulated. The number of orbits 
classified as PHZ analytically divided by the num- 
ber of orbits classified as PHZ numerically yields 
the percentage of PHZ, etc. Tab. [3l shows that 
the global correspondence between both methods 
is quite convincing, which can be considered a 
strong indicator that the behavior of the respec- 
tive HZs is modeled correctly. 

Also, Figs, m m & [71 indicate that most of the 
significant deviations between numerical and an- 
alytical results occur near the border of orbital 
instability, especially for high mass and small pe- 
riod ratios. This might be expected for AHZs 
and EHZs given the approximations involved in 
determining the analytic estimates. In the case 
of PHZs, however, semi-analytical results suggest 
caution in using simulation outcomes as reference 
values. 

8. Results 

Diverse trends in the behavior of the different 
types of HZs can be seen in Figs. [H & [71 While 
the AHZ is almost independent of the binary's 
eccentricity and coincides well with the KHZ by 
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KWR93 for distant stellar companions, the PHZ 
and EHZ shrink with higher binary eccentricities. 
The fact that the PHZ and EHZ contract around 
the center of the KHZ emphasizes the importance 
of the changes in the planet's eccentricity, as the 
secondary's radiation alone could only account for 
one-sided features towards the outer edge of the 
KHZ. The PHZs seem most affected by strong 
perturbations and shrink to almost half their size 
before the systems investigated become unstable. 
Interestingly, in close S-Type II binaries with low 
eccentricity the extent of the EHZs and AHZs can 
reach beyond the predicted values by KWR93. 
This can be seen in Fig. [3 where a zoom on the 
outer border of the KHZ in a G2-F0 configura- 
tion is shown. Here an extension towards the sec- 
ond star of about 0.1 AU seems possible for the 
system's AHZ if the binary's eccentricity remains 
below the system's stability limit of eb — 0.2. 

9. Conclusions 

In this work, the impact of the second star on 
Habitable Zones in S-Type binary star systems 
with different stellar constituents has been inves- 
tigated analytically as well as numerically. The 
radiative contribution of the secondary on a ter- 
restrial planet is negligible in all but S-Type II sys- 
tems, if orbital stability of the planet is required. 
The gravitational influence of the second star on 
the other hand perturbs in the planet's eccentric- 
ity, which in turn can lead to substantial changes 
in planetary insolation. Therefore the secondary 
has indeed to be taken into account when calcu- 
lating the extent of Habitable Zones. 

Our analytical estimates for planetary eccen- 
tricities in binary star systems introduced in Ap- 
pendix |B] are an exten sion to secular perturb ation 
theory as used in e.g. Heppenheimeil ( 19781 ). To- 



gether with methods presented in section |6] sug- 
gesting an analytic determination of Habitable 
Zones, they allow to paint a global picture of hab- 
itability in S-Type binary star systems without 
having to rely on time consuming numerical or- 
bit integrations. Our approach is quite flexible 
in the sense that different planetary atmospheric 
models and average stellar luminosities can be in- 
tegrated via adaption of the Seff values. Thus, 
the formulae presented in this article grant access 
to calculating Habitable Zones for a large set of 



possible binary-planet conflgurations. 

For the three stellar conflgurations investigated 
it could be shown that the Permanently Habitable 
Zone, i.e. the zone where the planet never ex- 
ceeds the classical insolation limits for habitability 
shrinks considerably with the binary's eccentricity. 
If one considers average insolation values only, the 
extent of the Average Habitable Zone coincides 
well with predictions by iKasting et al. I (Il993l) for 
wide binaries, whereas a significant extension to- 
wards the secondary is possible for close, eccentric 
binary systems. The overall correspondence be- 
tween numerical and analytical results presented 
is excellent, as 93-99% of all investigated orbits 
were classified identically. The computational ef- 
forts required to calculate the true extent of Per- 
manently Habitable Zones numerically, however, 
can be enormous and might in fact be prohibitive 
in some cases. In contrast, the analytical method 
presented offers immediate, reliable estimates. 

A more careful approach than the one proposed 
in this work is advisable when multiplanetary sys- 
tems, systems close to the stability limit or res- 
onant configurations are being investigated. In 
a next step we plan to extend our classification 
methods to mutually inclined systems. 
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A. Insolation Variance 

The insolation variance a planet receives on an S-Type orbit in a binary star system was defined in section 
16.31 as: 

a' = {Sl,)t ~ {Stot)l (Al) 
Considering the linearity of the expectation value operator, the term (Stot)'^ can be decomposed irf^: 

{Stotf = (^i)' + 2{Si){S2) + {S2? (A2) 

Averages for insolation values from both stars have been derived in section IH] already, and are therefore not 
repeated here. Instead, we will develop expressions for the first term on the right hand side of equation (jAip . 

{SD = (Sf) + 2{S^S2) + {SD (A3) 

Using equivalence radii r and R for the planet and the secondary respectively, which were introduced in 
section m we get: 

(Sf) = ^ ihdt - dM = ^ = --A-^ (A4) 



^^^^'^ ^2^1 '^R^ + r^-Rrcos{cl>)'^'^ = R^r^ - ^^^^ 

~ 27tJ, [R^+r^-Rrcos{cly)) ~ {r^ - R^^ ^ ' 

where S(t) is again the time dependent distance of the planet to its host star, and Up and Cp the planetary semi- 
major axis and eccentricity respectively. Given the circular nature of the orbital equivalence approximations 
we have applied, the results are bound to underestimate the true variances. Since the radiative contribution 
of the primary dominates the planetary insolation for S-Type I systems and is at least as improtant as the 
secondary's insolation for the S-Type II systems investigated, we will use stronger estimates for {Sf)t than 
given in relation (|A4[) : 

r2Tr 



Lin /l + (ep)cos(/)^' 



df 



As one can see, the difference between relations (jA4p and (jA7p is negligible for small injected planetary eccen- 
tricities, but its contribution becomes important if the injected eccentricities grow. Combining expressions 
l7|) , (|A5I) and (jA6p with the respective terms of {Stot)t produces the desired variance: 



a' 



2\2 



3 ( -1 + H4) Helf + {elf + ( 1 - %^ ^'^^ 



^We drop the subscript t on the averages, as there is no danger of misinterpretation. 
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B. MELximum Planetary Eccentricity 

The maximum possible eccentricity a terrestrial planet's orbit can acquire in an S-Type setup is composed 



of 



where e^^ denotes amplitude of short-period t erms and e''^'^ the planetary eccentricity's secular amplitude. 
Using expressions derived in iGeorgakarakosI (|2003l ) via the Laplace-Runge-Lenz vector, we can estimate 
the maximum short-period contributions for the planet's eccentricity by taking only terms including the 
dominant fre quencies into account . The amplitude of the secular part of the planet's eccentricity is given by 
2C/{B - A) (lGeorgakarakosll2003h resulting in the following expressions: 



a 



15 13 (4-|-lle2) 11 1 (l + eb)3 3 1 {I + eb)\& + llcb) 



64 (1 - e 
b a 3 



2)5/2 

2el 



4 (1 
2 



4^3 



(l-e2)9/2 



4X1/3 (1 -e2)i/2 



7X1/3(1 _e2)i/2 + ±x2/3(l-e^) 



(B2) 
(B3) 



The mass parameters a,/3 and 7 are defined as: 



m2 



rriinip 



Ml/3 



(mi + mpY'^M^/^ 



m2{mi + nipY/'^ 



nip being the planetary mass and mi,m2 the stellar masses of primary and secondary respectively. 

M = TOi + 1712 + rup is the total mass of the system. Finally, X denotes the secondary to planet period ratio: 



X 



mi 



M 



1/2 



3/2 
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Fig. 3. — toy left: Higliest possible binary eccentricity as a lunction of the secondary's semi-major axis, 
if a terrestrial planet was to remain dynamically stable on orbits with semi-major axes corresponding to 
the inner (runaway gr eenhouse, open symb ols) or outer (maximum greenhouse, full symbols) boundaries of 
the HZ as defined by iKasting et al.l ( 1993[ ). The curves represent quadratic fits of the FLI stability data 
presented in Tab. [2] Three different S-Type stellar configurations are shown: G2-F0 (A, A), G2-G2 (□,■) 
and G2-M0 (o,»). fop right: Minimum distance between planet and secondary (D) permitting planetary 
orbital stability in units of planetar y aphelion di stances from the primary (d = ap(l -I- e™'^^)). The planet's 
eccentricity was estimated following iGeorgakarak os (2005). bottom left: Here, the increase of the primary's 
effective insolation onto a terrestrial planet with injected eccentricity is compared to a single-star setup where 
the planet remains on a circular orbit. The planet is considered to be in periastron position (ap(l — e™"^)). 
bottom right: The secondaries' maximum radiative contributions to planetary insolation are presented - the 
planet is in apocenter position with regard to the primary (d), and the secondary is at pericenter. Once again, 
the results are normalized with regard to values which a single-star configuration with a terrestrial planet 
on a constant circular orbit would exhibit. One can see that the primary's radiative influence dominates in 
S-Type 1 systems, whereas for close S-Type II configurations the secondary's contribution is almost equally 
important. 
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Fig. 4. — Evolution ol insolation onto an Earth-like planet in a G2 - G2 binary star system. The oscilla- 
tions (left) are due to the injected changes in the planet's eccentricity (right) caused by the gravitational 
perturbations of the secondary. See the electronic edition of the Journal for a color version of this figure. 




Fig. 5. — Classification of HZs in a G2-F0 binary star system of S-Type II with a semi-major axis at = 10 AU. 
left: Zoom on the outer limit of the classical HZ (dashed line) close to the instability region (purple). The 
results were gained using analytic estimates presented in section El Black (blue online) denotes the PHZ, 
dark gray (green online) the EHZ, light gray (yellow online) the AHZ and white (red online) indicates 
that the planet is not habitable. The gray striped area (purple online) corresponds to the dynamically 
unstable region (IIW99), the striped extension shows the onset of dynamical chaos (PLD02). right: The 
numerical simulation results for the same configuration. The HZ limits extend beyond the values defined in 
Kasting et al. I (|l993l) . However, the white dashed-dotted line corresponds to the semi-analytic estimates of 



the PHZ using numerically derived values for e™°^. The semi-analytic results agree with the fully analytic 
estimates. This may indicate shortcomings of the entirely numerical approach to identify PHZ boundaries. 
The resolution of these calculations is Aap = 0.002 AU, and Aef, ~ 0.002. 
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Fig. 6. — Comparison of HZ classifications in three binary star systems with a secondary's semi-major 
axis of flfc = 50 AU. left: analytical estimates as discussed in sectionlH right: reference results gained 
via numerical integration, with a resolution of Aop = 0.01 AU, and Aeb — 0.01. The investigated stellar 
spectral configurations are: top: G2-F0, mid: G2-G2, bottom: G2-M0; The PHZ is represented in black 
(blue online), the EHZ is dark gray (green online) , light gray (yellow online) indicates the AHZ and white 
regions (red online) mean that the planet is outside of any defined HZ. The gray striped area (purple online) 
denotes dynamically unstable parameter regions (HW99), whereas the striped extension highlights the onset 
of dynamical chaos (PLD02), see section [3] The borders of the classical HZ as defined in KWR93 are 
represented by the vertical solid and dashed lines. 16 
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Fig. 7. — Same as Fig.[6]but for binary configurations with semi-major axis of af, — 10 AU . A decrease in the 
binaries' semi-major axes leads to more pronounced differences between analytic estimates and numerical 
simulation results. This can be expected from the approximations used to calculate PHZ, EHZ and AHZ, 
see section |6l Strong perturbations near the area of instability (gray striped - purple online) are modeled 
less accurately by the analytical estimates. The white dashed-dotted line corresponds to the semi-analytic 
estimates of the PHZ using numerically derived values for e™"^. 
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